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ABSTRACT 

The hermitian Wilson kernel used in the construction of the domain-wall 
and overlap Dirac operators has exceptionally small eigenvalues that make it 
expensive to reach high-quality chiral symmetry for domain-wall fermions, or 
high precision in the case of the overlap operator. An efficient way of suppress- 
ing such eigenmodes consists of including a positive power of the determinant 
of the Wilson kernel in the Boltzmann weight, but doing this also suppresses 
tunneling between topological sectors. Here we propose a modification of the 
Hybrid Monte-Carlo algorithm which aims to restore tunneling between topo- 
logical sectors by excluding the lowest eigenmodes of the Wilson kernel from the 
molecular-dynamics evolution, and correcting for this at the accept/reject step. 
We discuss the implications of this modification for the acceptance rate. 
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I. INTRODUCTION 



The domain-wall-fermion (DWF) [l|, 0] formulation of lattice QCD is based on a her- 
mitian, four-dimensional, super-critical Wilson-like kernel H\y = 7 5 Djy, which can be in- 
terpreted as governing propagation into an extra, fifth lattice dimension with coordinate 
s = l,...,/^. 1 For sufficiently large L 5 , right-handed and left-handed four- dimensional 
quark fields live on walls opposite of each other across this fifth dimension, and are rep- 
resented by wave functions which are bound to these walls. As long as L5 is finite, the 
right-handed and left-handed quarks are coupled through a "residual" mass [1, 0, @|, an 
effective chiral-symmetry breaking mass term originating from the overlap of their wave 
functions. 

In the free theory [l|, [6| or in perturbation theory jl, 0] these wave functions are expo- 
nentially bound to the walls, and therefore the residual mass is exponentially small. It is 
proportional to e~ aminL5 where the gap of the transfer-matrix hamiltonian, a m i n , is non-zero 
whenever the gap of the free Hy/ is. On realistic gauge field configurations, however, the 
super-critical H w has a spectrum reaching all the way down to zero (1, 0], with a mobil- 
ity edge that, by definition, separates a low-eigenvalue spectrum of localized modes from 
a high-eigenvalue spectrum of delocalized, or extended, modes. The near-zero modes are 
localized with a range of order the lattice spacing a, provided that the bare gauge coupling 
g and domain-wall height M are chosen to be away from the reg ion in the phase diagram 
where the mobility edge is less than order one in lattice units [10j |. 

An additional contribution to the residual mass is coming from the near-zero modes. It 
can be approximated by the integrated spectral density jlyl 5 d\ p(A), where p(A) is the 
ensemble- average spectral density of Hw (see Appendix C of Ref. |5|). If the near-zero 
spectral density depends only mildly on A, this contribution is roughly equal to p(0)/L 5 . 
Therefore, as L 5 is increased, eventually the slowly falling contribution of the near-zero 
modes becomes dominant. Much effort has been directed toward finding methods to suppress 
their presence, i.e., to make p(A) as small as possible for |A| < 1/L 5 . The same goal is relevant 
for overlap simulations, where suppression of near- zero modes is equally important in order 
to make it affordable to compute the sign function e(Hw) = Hw/\Hw\, which is at the heart 
of the construction of the overlap operator [111 , to high precision. 

A "surgical" method for suppressing the near-zero modes consists of modifying £>, the 



original domain- wall or overlap Boltzmann weight, to B det(H w ) [12l. Il3 ] . The first thing to 
note is that introducing the auxiliary determinant of H w into the Boltzmann weight does 
not change the universality class. Indeed, apart from a change in the lattice spacing, this 
modification does not affect the long-distance physics, because in itself the Wilson kernel 
H\y is a Dirac operator with a mass tuq of the order of the cutoff, \amo\ = 0(1). 

It has been demonstrated that, after re-adjusting the lattice spacing to the same value, 
the inclusion of detj Hw ) in the path integral results in a significant reduction in the near- 



zero mode density [12J, [13|. This is easily understood at an intuitive level because the 
Boltzmann weight now contains an extra factor of A 2 for each eigenmode of Hw, including 
in particular the low-lying, localized ones. Provided we stay away from the region where 



There is much freedom in choosing the super-critical kernel, and the detailed form is not important for 
this paper. For simplicity, we will assume the common choice that Dyy is the standard Wilson operator 
with bare mass m in the region — 2 < am < 0. The domain- wall height is M = |arao|. 
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the mobility edge is too small, the near-zero modes of H\y are associated with lattice-scale 
dislocations in the gauge field [§]. One may expect the inclusion of det(if^) to deplete the 
near-zero modes efficiently because this method suppresses selectively only those — relatively 
rare — dislocations that happen to support a near-zero mode. In other words, once we re- 
adjust the bare coupling such that the gauge field regains the same average "roughness" 
as before the insertion of det(if^), and hence the lattice spacing is unchanged, we find 
that only the "troublesome" dislocations leading to exceptionally small eigenmodes of H^y 
have been removed. The new spectral density behaves like p(X) ~ A 2 at small A, and the 
contribution to the residual mass is of order J^ul s ^ ~ V-^Ij t° be compared with a 
1/L 5 suppression for a spectral density that does not vanish at A = 0. 

In spite of their largely negative role discussed above, the near-zero modes also play 
what might be called a positive role. The gauge-field configuration space is divided into 
topological sectors whose boundaries are defined by configurations where det (Hw) =0. If 
we continuously follow a curve in the gauge-field space, it follows that the topological charge, 
as measured through the index of the overlap operator, changes by one unit precisely when 
one of the eigenvalues of Hw changes sign [141 ] . 

To date, the Hybrid Monte Carlo (HMC) algorithm [15| is the de facto algorithm for 
dynamical fermion simulations. The question is whether the trajectories generated by the 
HMC algorithm, or by one of its variants, will sample all topological sectors. Having fewer, 
even much fewer, near-zero eigenvalues would not constitute a problem all by itself, so long 
as these eigenvalues could move around and cross zero relatively freely. Unfortunately, the 
inclusion of det(H w ) not only reduces the near- zero spectral density significantly, but also 
completely suppresses the transitions between topological sectors. The reason is that, during 
the molecular-dynamics (MD) evolution phase of the HMC algorithm, the surfaces of co- 
dimension one where det(Hw) = constitute infinite-energy barriers. As a result, only 
one topological sector is sampled, typically the sector with zero topological charge. This 
constitutes a breakdown of ergodicity. In other words, the effect of adding in det(H w ) while 
using the standard HMC algorithm is to generate an ensemble according to the Boltzmann 
weight B det(H w )8(Q — q), where Q is the topological charge operator, and q the topological 
charge of the initial configuration. Thus the correct weight Bdet(H w ) is reproduced in just 
one topological sector, while the weight vanishes in all other sectors. 2 

Unlike 6*-vacua, fixed-topological-charge vacua do not cluster. This leads to enhanced 
finite- volume effects, which, already for single-particle masses, decrease only with an inverse 
power of the volume, and, moreover, increase with the inverse pion mass [16]. This must 
be compared with the usual exponential fall-off of finite-size effects in QCD with massive 
quarks. 

In this paper, we propose a modification of the HMC algorithm intended to restore 
ergodicity by allowing for zero crossings in the spectrum of H w even though det(if^) is part 
of the Boltzmann weight. The generated ensemble will then reproduce the correct Boltzmann 
weight Bdet(H w ) for configurations in all sectors. We will refer to this modification as the 
"Tunneling Hybrid Monte Carlo" or THMC algorithm. 

The key idea is that, during the MD evolution, the usual HMC fermion force associated 
with H w is replaced by the force derived from a modified operator of the form of H w +aQ^Q. 
The effect of the new term, aQ^Q, is to lift the (lowest few) near-zero eigenvalues of H w . 



2 We assume that there are no exactly massless quarks, and thus B > in all topological sectors. 
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This removes the infinite-energy barriers. As a result, the MD evolution can, and, we 
argue, will, visit all topological sectors. Using a trick reminiscent of renormalizat ion-group 
blocking, the original determinant det(if^) can be factorized as det(ifj^ + aQ^Q) times 
another term that corrects for the lifted near-zero modes. The original Boltzmann weight 
Bdet(Hw) is restored by incorporating the correcting factor in the accept/reject step that 
separates successive MD trajectories. While this procedure must lead to some reduction 
in acceptance, we will give arguments suggesting that that reduction may be affordable, 
given the physical merits of the generated ensemble. This paper does not as yet contain any 
numerical tests of the algorithms; these will have to be carried out in the future. 

This paper is organized as follows. In Sec. [Ill we set up the framework and explain the 
basic idea. In Sec. IIHI we discuss the conceptually simplest realization of the algorithm, in 
which the new term that lifts the near- zero eigenvalues, aQ^Q, is chosen deterministically. 
The resulting implementation of the algorithm may not be very practical because of its 
cost. A more affordable implementation of the algorithm, introduced in Sec. IIV[ is based 
on a stochastic choice of the new term. In Sec. |V] we discuss reasonable choices for the 
parameters of the algorithm, and we end with our conclusions in Sec. IVII Some technicalities 
are relegated to two appendices. 



II. THE BASIC IDEA 

We will consider a lattice gauge theory with partition function 

Z = j VU exp(-S g (U)) det{H^(U)) (2.1a) 

= J VUV^Vcf) exp(-S g (U) - S pf (U)) , (2.1b) 

S pf = <pH£(U)<j> , (2.1c) 

where S g is the (unspecified) gauge action, and once again Hw is the hermitian, super- 
critical Wilson operator. The QCD partition function should of course also include the 
physical domain-wall or overlap fermion determinant for each quark flavor. However, since 
we are concerned only with the effects of including the unphysical determinant det(H^), 
we have dropped the determinants representing the physical quarks. The continuum limit 
of the lattice theory (12.11) is thus a pure Yang-Mills theory. 

In Eq. (I2.1bj) the determinant of Hl^ has been rewritten as a pseudo-fermion partition 
function. If ergodicity were to hold, the HMC algorithm would generate an ensemble of 
configurations Ui with a probability measure 

V g {U) = Z- 1 exp(-S g (U)) det(flfUW)) > ( 2 - 2 ) 

by alternately updating the pseudo-fermion field <fi and updating the gauge field U. As 
explained in the Introduction, the standard HMC algorithm in fact fails to be ergodic because 
of the zero modes of H\y, and this is what we seek to amend with the THMC variant. 



Let us recall how the HMC algorithm normally works (for recent reviews, see Ref. |17j|). 
An updating cycle begins by generating a new pseudo-fermion field via a heat bath. A 
random complex vector £ is drawn from a gaussian ensemble with unit width, that is, with 
probability distribution P(£) oc exp(— followed by setting 

<P = H W (U)£ . (2.3) 
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Thus the probability distribution for <p is T > {4 > ) °^ ex P( — 4^ H^{U)(j)). A set of fictitious 
momenta is similarly drawn, corresponding to the free action S w = \'Y1 IX „ a 7T^ )At)0 - The 
gauge-field update then consists of two steps. First, the initial configuration of gauge field 
and conjugate momenta, {U, 7r}, is evolved along an MD trajectory by numerically solv- 
ing the classical Hamilton equations with some "guiding" hamiltonian TCmd- I n the most 
straightforward case, one takes as guiding hamiltonian 

Hmd = S^tt) + S g (U) + S pf {U) . (2.4) 

The gauge- field dependence of S p f comes from the operator H^{U) in Eq. (12. left . We have 
suppressed the dependence of S p f on the pseudo-fermion field <fi, which is kept fixed during 
the MD evolution. Hamilton's equations are numerically solved with some symmetric inte- 
grator with a non-zero step size St, for a number of steps Hmd — t md/St. This introduces 
"step-size" errors, which at the end of an MD trajectory are corrected for by a Metropolis 
test: the final configuration {W, tt'} is accepted, which means that W becomes the initial 
gauge-field configuration for the next cycle, or rejected, in which case the next cycle begins 
with the same initial gauge-field configuration U, with probability 

P accept = min{l,exp(W(7r, U) - H{ir / ,W))} , (2.5) 

where we have taken Ti = riuo- Assuming ergodicity, after thermalization the generated 
sets of configurations {Ui, 7Tj, <pi} follow the probability distribution 

V(U, tt, <f)) = Z- 1 exp(-H(U, tt, 0)) , (2.6) 

which, upon integrating over the pseudo-fermions and fictitious momenta, reduces to the 
probability distribution (12. 2p . 

In many applications, the randomness of the fictitious momenta and the pseudo-fermions 



at the start of each MD trajectory would be sufficient to make the algorithm ergodic [18 
However, in the case at hand this is hampered by the existence of zero modes of H\y- 
During the MD evolution Hmd(^,tt,0) plays the role of "energy." The pseudo-fermion 
action, which is part of TLmd(M, TT) 0)> becomes infinite whenever H W (U) has an exact zero 
mode. Because the initial energy is always finite, and since it is (approximately) conserved 
during the (discretized) MD evolution, (in practice) the MD evolution never passes through 
any configuration where det(H w (U)) = 0. 3 As explained in the introduction, this implies 
that only one topological sector is sampled. 

What one would like to do is to somehow remove the infinite-energy barriers separating 
topological sectors during the MD evolution. The basic philosophy of the HMC algorithm 
allows for this option: one is free to choose ri in Eq. (12.51) different from the guiding hamil- 
tonian Hmd- Eliminating the infinite-energy barriers necessarily means that H-md cannot 
be given by Eq. (12.41) if S p f is given by Eq. (12. left . But any "mistake" made by choosing a 
different guiding hamiltonian can be corrected for at the accept /reject step if the appropri- 
ate hamiltonian is used in Eq. (12.51) . It should of course be anticipated that the discrepancy 
between the two hamiltonians will lead to reduced acceptance, an issue that we will return 
to below. 



3 We note that it is not possible to regain tunneling by choosing St to be very large because this must lead 
to significant drop in acceptance and ultimately to integrator instabilities [17[ . 
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We thus wish to split the hamiltonian of the metropolis test into two parts: 



7~L — Ti-MD + 'Hcorr , (2.7) 

where contributions of near-zero modes are kept out of the guiding hamiltonian T~Cmd, and 
are accounted for by the correction TC CO rr, to be added in the accept /reject step. This can 
be done as follows. Borrowing from the way ultra-violet modes are integrated out in a 
Renormalization-Group (RG) block transformation for a bilinear action, one can show that 
the fermion determinant can be split as 

det(if^) = a~ n det(M) det(A) , (2.8a) 
M = H w + aQ ] Q , (2.8b) 
A' 1 = a' 1 + Q H w 2 Q ] . (2.8c) 

Here Hw and M. are N x N matrices, A is an n x n matrix, and Q is an n x N matrix, 
usually referred to as the blocking kernel in the RG context. For a proof of Eq. f)2.8p . see 
App.S 4 

Here we will choose the kernel Q such that Q^Q projects onto the lowest n eigenmodes 
of H w , either exactly or approximately. 5 Let \ipi), i = 1, • • • , n, be the eigenvectors of the 
first-order operator Hw with eigenvalues A« chosen such that {\f}™ =1 is the set of n smallest 
eigenvalues of H w . Now choose another set of n vectors \xi) such that 

= (ipi\Xj) ~ Sij , i,j = l,...,n, (2.9) 

and define to be the N xn matrix whose z-th column is the vector \xi)- Equations f!2.8bj) 
and (I2.8cl) then turn into 

n 

M = H w + aJ2\Xi)(x l \ , (2.10a) 

i=l 

(A-% = + (xi\Hw 2 \Xi) ■ (2.10b) 

Assuming Eq. (12.91) . the rightmost term in Eq. (j2.10aj) lifts the n lowest eigenvalues of H w , 
and the corresponding eigenvalues of M. are of order a. All the remaining eigenvalues of 
H w are more or less unaffected, if we choose a to be of order A^ +1 . 
Using Eq. fl2.8al) we may write 6 

det{H w ) = det(A) J VtfVtp exp(-S pf ) , (2.11a) 

S pf = ^M' 1 ^ . (2.11b) 

The THMC algorithm is based on choosing T~Cmd as in Eq. (12.41) with S p f given by 
Eq. (12. lib!) . In order to 

recover the desired Boltzmann weight (12.21) . we will include det(,4) 



4 Equation (|2.8p is true for an arbitrary choice of the kernel matrix Q and of the parameter a ^ 0. There 
is no restriction on n and N (curiously, it is valid even for n > N). Furthermore, Eq. (|2.8p holds in a 
limiting sense even if Hw has an exact zero mode. 

5 As discussed later in detail, we have in mind choosing n very small and, thus, n <C N . 

6 Following common practice to leave unspecified the overall normalization of a path-integral measure, we 
disregard from now on the normalization factor a~ n in Eq. (|2.8ap . 
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in the correction term exp(— H cor r), cf. Eqs. (I2.7P and (12. 5p . The definition of the THMC 
algorithm is not complete until a method for choosing the vectors is specified. There are 
two basic options: the \xi) can be chosen deterministically or stochastically, and this leads 
to the concrete implementations of THMC presented in Sec. II II I and Sec. IIV} respectively. 
In the latter case, exp(— Ti corr ) is not precisely equal to det(^4), as we will see. In both 
implementations the gauge field will be distributed according to the probability measure 



Using only the new form of the guiding hamiltonian, it is easy to see that all topological 
sectors will now be sampled. Any configuration at a crossing point between two adjacent 
topological sectors has an exact zero mode of H\y- But, so long as Eq. (I2.9p holds, Eq. (12.1 daft 
implies that the spectrum of M. has an 0(a) gap, even if Hw happened to have an exact 
zero mode. Thus, the excess MD energy associated with a crossing is of order 1/a. This 
is a finite number, because of our choice of a. This means that the initial MD energy 
has a considerable probability of being larger than the minimal MD energy needed for the 
crossing. As a result, the MD trajectory can pass through configurations where Hy/ has an 
exact zero mode. Typically, both at the start and end points of an MD trajectory the gauge- 
field configuration will be well inside some topological sector, so that neither det(A(U)) nor 
det(^4(W)) will be exceedingly small. This already suggests that there is no inherent reason 
that the acceptance rate should be intolerably small. We will discuss the issue of acceptance 
when introducing concrete implementations in the next two sections. 

III. DETERMINISTIC IMPLEMENTATION 

In this section we present the conceptually simplest implementation of THMC. It is 
defined by setting 

\Xi) = \A) , (3-1) 

at each time step during the MD evolution, where are the n lowest eigenmodes of 
on the updated gauge field. It follows that Zjj = 5ij in Eq. (I2.9p . and also that 

n 

M = H 2 w + aY,\^i)^i\ , (3.2a) 

i=l 

(A-% = (e^ + Ar 2 )^- . (3.2b) 
The representation of the partition function which is being simulated is thus 

Z = J VUV^V(j> det(^) exp(-S 9 - S pf ) , (3.3) 



with S p f given by Eq. f!2.11bj) . The guiding hamiltonian is given by Eq. (12.41) . and the 



accept /reject step is done with 

P accept = min |l,exp(W MD (7r,ZY) - H md {k' M')) det (l^)) } ' (3 ' 4) 

Comparing to Eqs. (12.51) and (12. 7ft . we see that this amounts to choosing 7i corr = 
— logdet(^4). We note that since Ai cannot be written in the form X^X with X local, 
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one updates the pseudo-fermions via <fi = r(M.(U))£ with some high-accuracy rational ap- 
proximation 7 r(x) ~ \fx. 

While conceptually the simplest, the deterministic implementation of THMC involves 
several costly ingredients. The n smallest eigenvalues of H w and the corresponding eigen- 
modes must be re-calculated after each MD step. For the calculation of the fermion force 
we have to monitor for level crossings during each MD time increment. 8 If no level crossings 
occurred, the calculation of the fermion force requires 

8x,fi{\A){A\) = {8xJA)){A\ + IV , *)(**,Ai(V'il) > ( 3 - 5a ) 
5 X J^) = -(H w - A,)" 1 (1 - I^X^I) (S Xtlt H w ) \A) , (3.5b) 
where the last equality follows from first-order perturbation theory, and 

8 x ,n = i \ U x>fl -^jf h.c.J . (3.5c) 

The numerical evaluation of & r JV'i) using Eq. (I3.5ap may be difficult in the vicinity of a 



level-crossing point. For a practical solution to this problem, see e.g. Ref. [19] . 

The situation is different if level crossing occurs between levels n and n + 1. In this 
case the set of n lowest eigenmodes changes discontinuously which, according to Eq. (I3.2al) . 
results in a discontinuity in the operator M.. The discontinuous change in S p f and, by 
Eq. (12.41) . in the MD energy, gives rise to a ^-function singularity in the fermion force. This 
singularity must be handled by a reflection/refraction step analogous to that discussed in 
Ref. [2l[. The discussion of Ref. [2l| addresses specifically the discontinuity of the overlap 
operator arising from a zero crossing in the spectrum of the corresponding kernel operator, 
here assumed to be H w . But the recipe of Ref. [2l[ is actually rather general, and can 
be easily adapted to the case at hand. This completes the definition of the deterministic 
implementation of THMC. 

We comment in passing that one might be tempted to consider the alternative of evolving 
continuously the n eigenmodes that were the lowest levels on the initial configuration. This, 
however, would constitute a breakdown of reversibility of the MD evolution whenever the 
set of n lowest eigenmodes on the final configuration is different from the n states that have 
evolved continuously from the initial set. 10 Thus, reflection/refraction steps are unavoidable. 

The cost of the deterministic implementation of THMC involves two components: first, 
the extra cost of each MD step, and second, the reduction in acceptance due to the omission 
of det(A) from the MD evolution. We discuss them in turn. 

In a dynamical overlap simulation, 11 in order to improve the accuracy of the evaluation of 
the sign function of Hw one usually computes n ov of the lowest eigenvalues of H w , alongside 



7 By construction, M. does not have any extremely small eigenvalues, and we expect that an accurate 
rational approximation of \J~M. is feasible with relatively few poles. 

8 We assume that for some n' > n, the n' lowest eigenmodes evolve slowly enough that one can track their 
evolution from time r to r + St. 

9 In an overlap simulation, it is known that the presence of roughly degenerate near-zero eigenvalues of 
Hw may lead to large fermion forces and to low acceptance rates. Using THMC may help indirectly, 
because the presence of det(H w ) will lead to much reduced density of low-lying eigenvalues. For a direct 
approach to solving this problem, see Ref. [2C|. 

10 We thank Urs Heller for this observation. 



li 



For a recent review, see Ref. [2/ 
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with their eigenfunctions and the variations defined in Eq. (13.51) . Thus, so long as we choose 
n in Eq. (12. lOj) to be less than or equal to n ov , there is no extra cost involved in using them 
for the deterministic THMC algorithm. 

The exception is the occurrence of level crossing between eigenvalues n and n + 1 of 
H^r which, as explained above, must be handled by a reflection/refraction step. This new 
reflection/refraction step represents an additional cost not present in an ordinary dynamical 
overlap simulation. We stress that the discontinuity encountered here results directly from 
the breakup of det(H^) in Eq. (12.8aj) . The discontinuity in det(.M) is matched by a dis- 
continuity in det(.A), such that det(Hyy) as a whole is continuous. Also, clearly, in this case 
there is no discontinuity in the corresponding overlap operator, because no eigenvalue of 
H\y has crossed zero. The cost of this new reflection/refraction step could be more tolerable 
than in the familiar overlap case, because the eigenfunctions of the crossing eigenvalues are 
both expected to be well localized, hence the cost of the necessary computation is likely to 
be independent of the volume. In addition, also the number n of eigenmodes kept out of the 
MD evolution is held fixed and small, independent of the volume (see also Sec. IVI|) . 

Omitting det(^4) from the MD evolution, or equivalently, performing the metropolis test 
with a hamiltonian 7i which is different from the guiding hamiltonian T~Cmd, will in general 
decrease the acceptance. In App. [B] we consider a crude model to estimate this effect. 
According to this model, when a single eigenmode is omitted from the MD evolution, the 
drop in acceptance is expected to be in the range of 1/3 to, at most, 1/2. When two 
eigenmodes are omitted, the drop in acceptance is in the range of roughly 50% to 65%. It 
is thus desirable to keep the number of eigenmodes relegated to det(*4) as small as possible, 
perhaps even n = 1. We return to this issue in Sec. [V] In order to avoid confusion, we stress 
that in the stochastic implementation discussed in the next section there is an additional 
source for a reduction in acceptance. 



IV. STOCHASTIC IMPLEMENTATION 

We now turn to a different version of THMC in which the vectors \xi) are kept fixed during 
the MD evolution. As we will see, a valid algorithm exists provided that the vectors are 
chosen stochastically. While the n lowest eigenmodes of H%y will have to be computed at the 
beginning and the end of each trajectory, this implementation avoids the costly operations 
of re-computing these eigenmodes and their gauge-field derivatives at every MD step. Also, 
no reflection/refraction steps are needed for the trivial reason that the vectors \x%) are not to 
be evolved. 12 The stochastic implementation is thus more suitable for domain- wall fermion 
simulations, where these calculations are normally not done. 

Before we proceed, we need to address a technical issue. Eigenfunctions of H w are 
determined only up to an arbitrary overall phase. In Sec. IIII[ the undetermined overall 
phase dropped out (see in particular Eq. (13. 2p ). and this ambiguity thus was of no concern. 
For the stochastic implementation introduced below, the phase ambiguity will have to be 
resolved. If is an eigenmode with some arbitrary overall phase, we may for example 



Of course, in an overlap simulation a reflection/refraction step is in principle still needed at a zero crossing 
of an eigenvalue of the kernel [2^ ]. 
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choose 

\i><) = ^7 m ■ (4-i) 

Here |1) stands for the vector with all component equal to one. 13 It is easily seen that is 
still a normalized eigenvector, and that it is unchanged if we multiply |^) by some arbitrary 
phase. 

We are now ready to introduce the stochastic implementation of THMC At the start of 
an MD trajectory we compute, as before, the lowest n eigenvectors \ipi(U)), which are now 
free of any phase ambiguity by virtue of Eq. (14. ip . We then draw n random vectors \r]i) 
from a gaussian ensemble with probability exp(— 'yfjf]), and set 

\ Xi ) = + h) • (4-2) 

Equivalently, the vectors \xi) can be thought of as n new complex bosonic fields (each 
carrying the same set of indices as a Wilson-fermion field) with their own Boltzmann weight 
exp(— Sker), with kernel action 

n 

Sker = 7 £ (Xi - MUM - V>i(«)> • (4.3) 

i=l 

With the proper (gauge-field-independent) normalization we then have 

J V X ] V X eMSker) = 1 , (4.4) 

which we may thus insert into the original partition function Eq. (12.1 1) , obtaining (with 
Eqs. flXgaD and (gXED ) 

Z = J VUV$V<\)V^Vx exp(-S g - S pf - S ker ) det(A) . (4.5) 

The relation between the representation (14.51) of the partition function and the stochastic 
implementation of THMC parallels the relation between the representation (I3.3P and the 
deterministic implementation. At the beginning of each trajectory a new set of pseudo- 
fermions 0, vectors \Xi), and fictitious momenta ir are drawn, each from the corresponding 
heat bath. The MD evolution is then carried out with the same guiding hamiltonian as 
before, cf. Eq. (12. 4p . Finally, the metropolis test (12. 5p is done with the hamiltonian 

H = Hud + S ker - log det(^) . (4.6) 

Note that both Hmd and det(^4) are now functions of the (fixed) vectors \xi), but not 
directly of any eigenmode of H^y. The initial and the final lowest n eigenmodes, \ipi(U)) and 
\ipi(U')), are only required for the evaluation of Sker- 

The stochastic implementation of THMC generates the correct equilibrium distribution 
(12. 2p . for the same reasons that the original HMC algorithm works [15). The algorithmic 



13 Other choices for the reference vector can be made in case that the inner product would happen 

to be numerically unstable. 
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role of the new stochastic degrees of freedom, \xi), parallels the role of the pseudo-fermions 
in the standard HMC algorithm (or, for that matter, in the deterministic implementation 



of THMC). In particular, satisfying the detailed balance condition [15| requires that the 
classical MD evolution be reversible. This means that if {U, tt} is evolved into {W, tt'} then 
{W, — tt'} is evolved into {U, —tt}, when all other degrees of freedom are held fixed. This now 
includes both and \xi)- Just like in ordinary HMC, using a symmetric integrator ensures 
the reversibility of the MD evolution in the stochastic implementation of THMC. 

According to Eq. (14. 2 p the vectors \x%) drawn from the heat bath depend on the initial 
gauge field. One might infer from this equation the erroneous conclusion that "reversibility 
is not satisfied." If this were true, reversibility would not be satisfied in the original HMC 
algorithm in the first place, because according to Eq. (I2.3P the pseudo-fermion field drawn 
from the heat bath depends on the initial gauge field as well. What drives this confusion 
is that, as a whole, the HMC algorithm is not reversible; it generates a Markov chain that 
converges to a specified probability distribution irrespective of initial conditions, hence irre- 



versibly [la, [17|, [18| . The same goes for both versions of THMC. Reversibility is only required 
for the classical MD evolution, where, as stated above, it amounts to the interchangeability 
of {U, tt} and {W, — tt'} for any given, fixed set of values of the remaining degrees of freedom. 
Thus, reversibility of the MD evolution is guaranteed by the use of a symmetric integrator 
in the two versions of THMC, just like in standard HMC. 



V. CHOICE OF PARAMETERS 

The THMC algorithm has not been tested yet, and at present it is not known how well it 
works in practice. Assuming some particular choice of the physical-sector action (including 
the gauge action and some action (domain-wall or overlap) for the physical quarks) there 
are five parameters that determine the performance of the stochastic implementation of the 
algorithm: n, a, 7, r MD and 5t. In Sec. |TT]we already argued that the "blocking" parameter 
a should be chosen equal to a typical small eigenvalue of H^. In this section we discuss the 
remaining parameters of the stochastic implementation. 

We begin with the parameter 7 occurring in Sk er , and explain why this parameter should 
be chosen to be of order one. We will assume that the initial and final low eigenmodes are 
all localized and, thus, that their components are of order one within their region of support. 
The norm 14 of the random vector \rj) is of order 1/1/7. We immediately see that one cannot 
choose 7 too small, if we want \xi) to be a reasonable approximation of \ipi{U)). If not, 
there would be no reason for the n lowest eigenvalues of to be lifted by the addition of 
aEi\Xi)(Xi\, cf. Eq. <^JQ^. 

The parameter 7 cannot be chosen arbitrarily large either. In fact, in the limit 7 — ► 00 
acceptance will vanish. For 7 ^> 1 the vectors \xi) drawn from the heat bath at the beginning 
of the trajectory will satisfy \xi) = \ipi{U)) + 0(1/-^). At the end of the trajectory we will 
thus have 

7 ( Xi - MW)\xi - = 7 MM - Wu')MU) - MW)) , (5-i) 

where we have neglected subleading terms in I/7. The accumulated discrepancy 



Because I77) is not localized, the magnitude of its individual components is of order 1/y/jV, where V is 
the number of lattice sites. 
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(ipi(U) — ipi(W)\il)i(U) — ipi(U')) is a function of various simulation parameters (in partic- 
ular the total MD time tmd), but it is independent of 7. On the final configuration we 
would thus have Sk er — > 00 for 7 — > 00, implying that every trial configuration would be 
rejected. Combining both requirements, we must therefore choose 7 of order one. 

There are three more parameters to discuss: the number of "lifted" eigenvalues n, the 
total MD time per trajectory tmd, and the MD time increment Sr. While existing numerical 
simulations of the super-critical det(iif^) using the ordinary HMC algorithm are confined 
to the sector with topological charge zero, the results can guide us in the choice of the 
remaining parameters. Our conclusion is that we advocate first trying THMC with values 



of tmd and 5t similar to those used in Refs. |12|, [13|, [23j, and with a very small value of n, 
perhaps even n — 1. 

The first interesting lesson is that, even at rather large lattice spacings (with a -1 as low 
as 1.5 GeV), near- zero eigenvalues of Hy/ were efficiently suppressed by the inclusion of 
det(ifj^) in the Boltzmann weight. What this means is that, even if we were to choose 
an unrealistically large value for n {e.g., n = 100), topology change during a given MD 
trajectory should still take place mostly through the zero crossing of just one, or very few, 
eigenvalues. In principle, then, THMC will restore ergodicity even if we allow for just one 
crossing per trajectory, which, in turn, is made possible by choosing any non-zero value 
for n. Because multiple crossings per trajectory should be much less frequent, a plausible 
strategy is to make no special effort to allow for them. 

The next question is about the effect of the number n of eigenmodes kept out of the MD 
evolution on the acceptance rate. As explained in Sec. [Til when using the THMC algorithm 
the acceptance rate will be less than one even in the limit 5r — > where HmdKjW) = 
rihini 71 ">^0- I n the deterministic implementation of Sec. 1111} acceptance is controlled by the 
ratio det(^4(W))/det(v4(W)) in the limit 5r — > 0. Based on the model of App. [B] we have 
concluded that the acceptance rate is likely to go down with increasing n, and thus it is 
desirable to keep n as small as possible. 

In the stochastic implementation, the metropolis test involves the additional factor 
SSker = SkeriU') — Skerip() (c/. Eq. (I4.6P ). Let us discuss the role of SSk er at a qualitative 
level, considering first the case n — 1. Because 7 = 0(1), we clearly have Sk er (p() =0(1) 
at the beginning of the trajectory. During the MD evolution the shape and location of the 
(localized) eigenmode |?/>i) changes by an unknown amount, and at the end of the trajectory 
Sker(U') will assume another 0(1) value. If we would now allow for larger values of n, the 
size of the interval containing the most probable values of 5Sker will increase. Positive values 
of SSker, which one expects to occur a fraction of the time, will thus tend to further decrease 
the acceptance rate with increasing n. We conclude that one should indeed try to keep n as 
small as possible. 

This brings us to the trade-off between tunneling and acceptance. While it remains to 
be tested, we think that it is not implausible that THMC with values as small as n = 1 or 



arc 



n = 2 will produce an appreciable tunneling rate. Many of the results of Refs. [12|, [13 
relevant for this question, but particularly illuminating is Fig. 2 of Ref. [23[ which reports the 
results of an experimental two-flavor simulation. We consider the left panel, generated with 
det(iffy) in the Boltzmann weight. 15 This plot shows the lowest few eigenvalues of Hw as a 
function of MD time. The MD time interval corresponding to one trajectory was tmd = 0.5, 



15 We note that Refs. 



13|,|23[ actually used det(i?^ / )/det(if^ / + /1 2 ) with an = O(l), in order to reduce the 
contribution of the large eigenvalues of H^y to the fermion force. 
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with tmd/^t = 0(20). We learn from the plot that, frequently, Hw has a small eigenvalue 
that changes by nearly 100% over a one-trajectory interval. This is an encouraging result, 
suggesting that if one would switch from HMC to THMC, there is a good change that zero 
crossings will take place. 

Two worries immediately come to mind. Consider those trajectories where one eigenvalue 
has changed by a relatively large amount and became exceptionally small toward the end 
of the trajectory. As already noted, such trajectories are relatively frequent in the data. 
However, it is often true that that particular eigenvalue is not the smallest, nor even the 
next-smallest (in absolute value), at the beginning of the trajectory. If we would now re-run 
the trajectory from the same starting point with THMC, but use only n = 1 or n = 2 and 
no higher value, our selection of the initial lowest eigenvalue(s) would miss that particular 
eigenvalue, which is the best candidate for crossing. In such a case, THMC would not do its 
job, because the \xi) would likely not have much overlap with the eigenvector corresponding 
to this eigenvalue, cf. Eq. f)2.9p . Nevertheless, for some fraction of the configurations the 
best-candidate eigenvalue did actually start as the smallest one. Thus, even n = 1 should 
produce a non-vanishing tunneling rate. 

A related worry is that when an eigenvalue became exceptionally small at the end of 
the trajectory, often the final configuration was rejected. 16 Once again, part of the answer 
is that some configurations were accepted even though they had one exceptionally small 
eigenvalue. If we would switch to THMC for the next trajectory, that eigenvalue would now 
be selected even with n = 1, and would have an appreciable chance to cross zero. 

However, more importantly, with ordinary HMC, any eigenvalue must eventually be de- 
flected away from zero because of the unboundedly growing fermion force. In contrast, with 
THMC any large repelling force originating from one of the n eigenvalues which were the 
lowest ones at the beginning of the trajectory is eliminated. An eigenvalue motion that, 
with standard HMC, had been slowed down by the unboundedly growing repelling force, 
will speed up under THMC. Therefore, we expect that many of the eigenvalue flows that 
ended up very close to zero (but with, invariably, no crossing) when HMC was used, would 
move more rapidly when we switch to THMC. In fact, they may well cross zero during the 
same MD time interval, and end up with the opposite sign and well away from being ex- 
ceptionally small by the end of the trajectory, thus increasing considerably the acceptance 
probability of the configuration at the end of the trajectory. 

Finally, we consider how the location of an initially low-lying mode evolves. If the eigen- 
mode would "drift" appreciably away from its original location, this would give rise to in- 
creased Sker on the final configuration, and once again decrease the acceptance rate. Worse, 
the overlap between the and this mode would rapidly decrease, thus possibly decreasing 
the smallest eigenvalues of A4, and hence re-introducing barriers in the MD evolution. This, 
however, is unlikely to happen for the following reason. 

Having an exceptionally low eigenvalue requires two things: values for some plaquettes 
substantially different from one, as well as a certain "fine tuning" of the link variables relative 
to each other. This picture, put forward in Ref. [§], is corroborated by the fact that large 
plaquette values are far more frequent compared to exceptionally low eigenvalues. Now, 
because the momenta conjugate to each link variable are basically uncorrelated with each 



This can be seen from the discontinuities in the spectral flows at integer values of the horizontal axis 
("traj"): when a configuration was rejected the next trajectory was started with the same initial config- 
uration as the previous one. 
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other, the likely outcome of the MD evolution is that any finely-tuned relation between 
local plaquettes will be quickly destroyed. Once again this expectation is supported by 
Fig. 2 of Ref . [23( , which shows that exceptionally- low eigenvalues tend to evolve much more 
rapidly compared to the "bulk" of (not-so-small) eigenvalues. This supports the conjecture 
that indeed the necessary fine-tuned dislocations only exist for a relatively short MD time, 
in agreement with the general observation that dislocations supporting exceptionally low 
eigenvalues are rare. In contrast, a "drift" of the location of the eigenmode (in which the 
eigenvalue is kept very small) would require an orchestrated motion of the dislocation with 
almost no change in its shape, something which is very unlikely to happen. 

The picture that emerges, then, is that of a localized eigenvector that remains more or less 
in place and localized, with a rapidly changing eigenvalue and, presumably, correspondingly 
rapid changes in its detailed shape. This implies that (ipi\Xi) remains of order one, and thus 
no MD barrier is re-generated. We emphasize, as above, that this picture only needs to hold 
for a sizable fraction of the low-lying eigenmodes that exist at the start of the MD evolution 
trajectories. 



VI. CONCLUSION 

In this paper we proposed a new variant of the HMC algorithm. The aim of the tunneling 
HMC algorithm is to allow for eigenvalue crossings through zero even in situations where 
the determinant containing such eigenvalues is included in the Boltzmann weight. The 
application we considered explicitly is lattice QCD with domain-wall or overlap fermions and 
with the auxiliary determinant of H^. The THMC algorithm may find other applications 
as well. 17 

A major issue is the competition between tunneling rate and acceptance. We have dis- 
cussed in some detail why we believe that appreciable tunneling is not incompatible with 
reasonably high acceptance. Setting aside the details, the underlying reason why this can 
be true is that the difference between the guiding and metropolis hamiltonians, TC — Hmd, 
is not an extensive quantity. For any lattice volume, in order to allow for tunneling between 
topological sectors, we only need to keep a fixed, and very small, number of degrees of free- 
dom outside of the guiding hamiltonian. Once ergodicity is restored via THMC, the correct 
Boltzmann weight will automatically be reproduced, in all topological sectors. Therefore, the 
cost of THMC does not grow with volume. We stress that our semi-quantitative conclusions 
will have to be tested numerically. 

We discussed two concrete implementations of THMC, a deterministic and a stochastic 
one. We believe that the stochastic implementation of THMC is likely to be of most prac- 
tical interest (especially for domain- wall simulations), because of its smaller computational 
overhead. For this to be the case, any further decrease in acceptance coming from the extra 
element present in the stochastic implementation, Sk er , should not be too large. 

In comparison to standard HMC, the extra cost of THMC comes partly from its more 
elaborate structure. For example, it takes more operations to perform a single multiplication 
of a vector by A4, as compared to H^- However, it is likely that the major extra cost of 
THMC will come from reduced acceptance: if acceptance is a factor of two smaller, this 



17 



For example, we envisage using THMC as an ingredient in the simulation of the ghost sector of equivari- 



antly gauge-fixed Yang-Mills lattice theories 24 1 . 
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immediately translates into a factor of two extra cost. Any extra cost of THMC will have to 
be weighed against the physical merits of the generated ensemble. For example, the fixed- 
topology simulations of JLQCD [l3| give rise to enhanced finite-volume effects proportional 
to l/(m£V) 0. Ultimately, the only way to keep such enhanced finite-size effects under 
systematic control is to simulate bigger lattice volumes, something that carries its own price 
tag. Clearly, it is the cost for "equal-quality" physics output which must be compared. For 
domain-wall fermion simulations, the extra cost of THMC must be balanced against the 
improvement of chiral symmetry. 

Are there possible alternatives to THMC? The answer is yes. A simple approach is to 
replace det(if^) by det(H w + e 2 ) as the auxiliary determinant, where presumably e 2 would 
be chosen comparable to a typical small eigenvalue of H^, much like the parameter a 
introduced in Eq. (12.81) . Because there are no exact zero modes any more, one can simulate 
det(Hyy + e 2 ) using ordinary HMC. 18 In comparison with THMC, this avoids altogether any 
drop in acceptance associated with a non-zero "H — 'Hmd- At the same time, if e is chosen 
small enough, the low-lying eigenvalues of Hyy would still be suppressed. An additional 
benefit of simulating det(H w + e 2 ) is that no new code needs to be written, and there are 
fewer new parameters that one has to experiment with for optimal performance. 

In spite of the obvious advantages of this alternative, we believe that it is an open question 
which solution is better. When using det(H w + e 2 ) one will have to compromise the value of 
e between two conflicting requirements: the basic need to suppress the low-lying eigenvalues, 
and the need not to suppress them too much so as to allow for tunneling. In contrast, the 
advantage of THMC is that it solves one problem at a time: the reduction of an abundance 
of low-lying eigenvalues comes from det(H w ), while the restoration of tunneling comes from 
the modification of the algorithm. 

While in this paper we have assumed that the auxiliary determinant is det(H w ), the 
THMC algorithm evidently generalizes to other even powers as this determinant. If the 
auxiliary determinant is det(H^), the near- zero spectral density will behave like p(A) ~ X 2k . 
In a domain-wall simulation, the residual mass will thus scale like 1/L 2fc+1 . This constitutes 
a bigger pay-back for any increase of L5 in the region where the residual mass is dominated 
by the near-zero modes. 

The choice of the JLQCD collaboration to run fixed-topology dynamical overlap simu- 
lations stems from several technical advantages arising from the introduction of det(if^). 
The (much!) smaller near- zero spectral density leads to a smaller range for l/\Hw\, and an 
affordable cost for an adequate approximation of the sign function Hw/\Hw\- The costly 
reflection/refraction step is not needed by design, because no eigenvalue of Hw ever reaches 
zero during the MD evolution. With THMC, the same reduction in the near-zero spectral 
density is achieved, with all the ensuing advantages. Of course, reflection/refraction steps 
cannot be avoided if we want topology to change. 19 But the number of attempted topology 
changes will remain small, and will not grow with the volume, because it is controlled by 
the number n of eigenvalues which are relegated to det(.4). As we have argued above, this 
number can be chosen to be very small. The cost of reflection/refraction steps should thus 
be tolerable. 

Some of the difficulties in achieving a satisfactory topology-changing rate are inherent. 
While too-many small eigenvalues of Hw is the problem at relatively large lattice spacing, 



We learned about this alternative approach from Norman Christ. 
See, however, Sec. 3.3 of Ref. [22j |. 
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when the lattice spacing gets smaller eventually the situation changes, because the near- zero 
mode density decreases rapidly with the bare coupling. This means that tunneling becomes 
rare, even without the auxiliary determinant in the Boltzmann weight. For other difficulties 
associated with tunneling from the Q = sector to sectors with a higher topological charge, 



see e.g. Ref. 22]. It is an open question whether THMC can help in creating a bigger interval 



of lattice-spacing values where the tunneling rate is satisfactory. 
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APPENDIX A: PROOF OF EQ. GZM 

We represent det(-ffjy) as a Grassmann integral over ip and if), multiply the integrand by 
unity written as a new Grassmann integral over \l/ and and perform the integration first 
over if> and if) and then over \1/ and Keeping all normalization factors the result is 

Z = det(H^) = J VipV^p exp(-^H^iP) (Al) 

/n 
Thf;Thp JJ d%d% exp {-^H 2 w if) - a(* - - QVO) 

/n 
Y[dVid$i exp(-C4#) 

i=l 

= a~ n det(Af) det(^) , 
where Ai is given by Eq. (12.8b|) and 

A = a- a 2 QM~ l Q ] . (A2) 
We next consider the two-point function 

/n 
Vif)V^]Jd%d% (A3) 



x exp (-if)Hlrif) - a(tt - ^Q f )(^ - Qif))) *< 
ThpVlji exp(-^V) (a -1 ^ + (Qip)i 



7-1 



Performing the integral over if) and if) we arrive at Eq. (I2.8cj) with the identification A = A. 
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APPENDIX B: ACCEPTANCE RATE IN THE DETERMINISTIC IMPLEMEN- 
TATION 



In this appendix we discuss a simple model for the decrease in acceptance resulting from 
the omission of det(^4) from the MD evolution, within the deterministic implementation of 
THMC We will first consider the case n = 1, namely, a single eigenmode is omitted from 
the MD evolution. We begin by writing down a general expression for the acceptance rate. 
At the beginning of a trajectory, x = det(A(U)) follows a probability distribution P(x). 
At the end of the trajectory, y = det(^4(W)) follows another distribution P(y\x). This 
distribution is conditional on x, because the initial and final gauge field configurations are 
not uncorrelated. Ignoring step-size errors, i.e., assuming that 71md(j{,U) = TCmd(^' M')i 
Eq. (13.41) gives an acceptance rate 



Pacce P t= / dx P{x) \ dy P{y\x) min{l, y/x} , (Bl) 
Jo Jo 

where have we rescaled the possible values of det(^4) to a unit interval. 

In order to estimate P aC cept we need to make several assumptions. The initial configura- 
tion U is an equilibrated configuration of the ensemble generated by the Boltzmann weight 
exp(— Sg{U)) rii^i' c f- Eq. f)2.ip . We assume that the eigenvalue Ao with smallest abso- 
lute value is localized, and practically uncorrelated with the remaining eigenvalues. This 
suggests a probability distribution for this eigenvalue which is proportional to Ag itself, for 
Ag <C 1. As for the support, we postulate a sharp-cutoff model, namely, we assume that 
the probability distribution for the smallest eigenvalue vanishes outside the interval [0, Cot]; 
inside this interval, it increases with Ag, starting from zero at Ag = 0. Physically, Ca is 
the average of the second-smallest eigenvalue (squared). Our choice of a (see Sec. [TT]) en- 
sures that C = 0(1). Treating C as a free parameter thus allows us to get an idea on how 
sensitively the drop in acceptance depends on the details of the theory. 

What enters Eq. (IBlj) is not directly the probability distribution for Ag but, rather, that 
for the related quantity x = det(A(U)). Treating x as the independent variable, we specify 
the probability distribution by postulating P(x) oc Xq(x) = xa/(a — x) over its support, 
where the last equality follows from Eq. (13.2bj) . Similarly using Eq. (13.2bl) to obtain the 



x-range that corresponds to the previously selected range of Ag, followed by a rescaling to 
the unit interval, we obtain the probability distribution 20 

P ^ = N C + l-Cx ' ° " X ~ 1 ' (B2) 

where 

^£±i,og(C + l)-I. (B3) 

We next turn to the conditional probability P(y\x). The small eigenvalues of on a 
typical configuration are expected to be well below its mobility edge, and thus to be localized 
with a range of order a (Tol . 25| . These eigenvalues are only sensitive to the details of the 



gauge field near the corresponding localized modes, and these gauge fields fluctuate on a 
short (MD) time scale. It is therefore not unreasonable to expect that, for typical values 



See below for more discussion of this crude approximation. 
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of the MD time interval tmd for one complete trajectory, the variable y actually becomes 
decorrelated from x, P(y\x) = P{y). Moreover, the THMC algorithm is designed to remove 
the lowest eigenvalue from the MD force. This suggests that the probability distribution 
P{y) is rather flat. For the sake of simplicity we thus take P(y) to be constant, and have the 
same support as P{x), which amounts to choosing P(y) = 1 for < y < 1 (after a similar 
rescaling). Finally, plugging this into Eq. (1B1|) and using Eq. (1B2I) we obtain 

r accept 2 ' ^(J ' 1° V 

The limits C — > and C — > oo give rise to P accep t = 2/3 and P aC cept = 1/2, respectively. 

Let us now generalize this model to the case that two eigenmodes are omitted from the 
MD evolution. We write det (A(U)) = x\x 2 , where x\ and x 2 are related to the corresponding 
eigenvalues in the same way that x was previously related to the smallest eigenvalue. These 
modes are expected to be localized, and uncorrelated, at the beginning of the trajectory. 
We may thus assume the same probability distribution as before, separately for x\ and for 
x 2 . 21 Similarly, writing det (A(W)) = y±y 2 , we will assume the same flat distribution for y 1 
and 2/2 at the end of the trajectory prior to the metropolis test. The acceptance rate is then 
given by 

Paccept = dx x dx 2 P{x x )P{x 2 ) \ dy 1 dy 2 min{l , {yxy 2 ) / \xix 2 )} , (B5) 
Jo Jo 

where P(x) is given by Eq. flB2j) . Varying C in the range of 0.01 to 100 we now find that 
Paccept varies between roughly 50% and 35% respectively. 

Returning to the case that just one eigenmode is omitted from the MD evolution, our 
model, the probability distribution P(x) of Eq. (IB2I) . neglects most factors in the Boltzmann 
weight. We have ignored the pure-gauge action, as well as all the remaining eigenvalues of 
H w . In addition, we neglected the jacobian that arises if we treat x = det(A(U)) as an 
independent integration variable. The motivation for neglecting all these factors is twofold. 
First, clearly, it would be difficult to incorporate all factors accurately. Also, while all the 
neglected factors are unlikely to be constant for the range where our chosen P(x) is non-zero, 
they are also not expected to vanish. The most crucial feature of the exact, yet unknown, 
probability distribution is that it vanishes with Aq. This essential feature is captured by 
our model. Of course, many other models with the same behavior for Aq 1 could be 
conceived. The free parameter C reflects to some extent the arbitrariness of the model. The 
fact that rather large changes in the value of C do not lead to a dramatic change in the 
resulting acceptance is, we believe, an encouraging sign that the acceptance rate will indeed 
be roughly in the range predicted by our model. Of course, ultimately the only way to verify 
these expectations is through numerical tests of the THMC algorithm. 
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